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ABSTRACT 

Hot, underdense bubbles powered by active galactic nuclei (AGN) are likely to play a key role in 
halting catastrophic cooling in the centers of cool-core galaxy clusters. We present three-dimensional 
simulations that capture the evolution of such bubbles, using an adaptive-mesh hydrodynamic code, 
FLASH3, to which we have added a subgrid model of turbulence and mixing. While pure-hydro 
simulations indicate that AGN bubbles are disrupted into resolution-dependent pockets of underdense 
gas, proper modeling of subgrid turbulence indicates that this a poor approximation to a turbulent 
cascade that continues far beyond the resolution limit. Instead, Ray leigh- Taylor instabilities act to 
effectively mix the heated region with its surroundings, while at the same time preserving it as a 
coherent structure, consistent with observations. Thus bubbles are transformed into hot clouds of 
mixed material as they move outwards in the hydrostatic intracluster medium (ICM), much as large 
airbursts lead to a distinctive "mushroom cloud" structure as they rise in the hydrostatic atmosphere 
of Earth. Properly capturing the evolution of such clouds has important implications for many ICM 
properties. In particular, it significantly changes the impact of AGN on the distribution of entropy 
and metals in cool-core clusters such as Perseus. 

Subject headings: hydrodynamics - cooling flows- Xrays: galaxies: clusters 



1. INTRODUCTION 

The X-ray and abundance profiles of the hot, diffuse 
medium in galaxy clusters are observed to be bimodal 
(Fukazawa et al. 2000; Matsushita et al. 2002; Schmidt 
et al. 2002; Churazov et al. 2003; De Grandi et al. 2004). 
Strong intracluster medium (ICM) abundance gradients 
are associated with cool-core clusters with a peak in their 
central X-ray surface brightness distributions, and nearly 
uniform abundance profiles are associated with non cool- 
core clusters. Furthermore, these metallicity distribu- 
tions are much broader than the associated galaxy light, 
indicating that significant mixing has occurred (e.g. Chu- 
razov et al. 2003; Chandran 2005; David & Nulsen 2008). 

At the same time, the nature of cool-core clusters re- 
mains uncertain. Although strong X-ray emission indi- 
cates that the central gas is cooling rapidly, the deficit 
of star formation and < 1 keV gas (e.g. Fabian 1994; 
Tamura et al. 2001) means that radiative cooling must 
be balanced by an unknown energy source. Currently, 
the most successful model for achieving this balance re- 
lies on heating from a central AGN, yet the details of this 
process are poorly understood (e.g. Loken et al. 1995; 
Bruggen & Kaiser 2002; Reynolds et al. 2002; Brighenti 
& Mathews 2006; Brunetti & Lazarian 2007). 

While AGN are observed to drive large bubbles filled 
with relativistic particles (e.g. Boehringer et al. 1993; 
Carilli et al. 1994; McNamara et al. 2000; Blanton et al. 
2001; Finoguenov & Jones 2001; Nulsen et al. 2005), the 
synchrotron radiation emitted by the electrons in these 
regions fades after about 10 8 years, becoming extremely 
difficult to detect. Moreover, the corresponding depres- 
sions in the X-ray surface brightness are only visible near 
the center of the cluster, where the contrast is large. 
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Thus, it is unclear how far these structures rise in the 
cluster. Furthermore, AGN have been observed to in- 
duce shocks and/or sonic motions in the ICM that are 
believed to eventually dissipate their energy into this gas 
(Fabian et al. 2003; Kraft et al. 2003; Ruszkowski et al. 
2004; McNamara et al. 2005), although the impact of the 
resulting heating is difficult to quantify observationally. 

The presence of AGN-heated cavities in galaxy clus- 
ters has raised a number of questions. These buoyant 
bubbles are unstable to the Rayleigh- Taylor (RT) insta- 
bility, which occurs whenever a fluid is accelerated or 
supported against gravity by a fluid of lower density. In 
any ideal hydrodynamic model, the cavities must be in- 
flated supersonically or else they would be destroyed by 
RT instabilities faster than they are inflated (Reynolds 
et al. 2005). Curiously, the strong and hot ICM shocks 
that are expected around the active cavities are absent, 
and, instead, many cavities are surrounded by shells of 
gas that is cooler than the ambient ICM (Nulsen et al. 
2002). Secondly, these cavities appear to be intact even 
after inferred ages of several 10 8 yrs, as for example the 
outer cavities in Perseus (Nulsen et al. 2002). However, 
hydrodynamic simulations fail to reproduce the observed 
morphology as the RT and other instabilities shred the 
bubbles in a relatively short time (Robinson et al. 2004; 
Reynolds et al. 2005, Bruggen et al. 2005a), although 
this time can be extended somewhat by a more detailed 
treatment of bubble inflation (Pizzolato & Soker 2006; 
Sternberg et al. 2008; Sternberg & Soker 2008a,b). 

A consequence of the evolution of such cavities is the 
production of turbulence, which is likely to be pervasive 
in the ICM and play a crucial role in the structure of 
cool-core clusters (e.g. Schuecker et al. 2004). Turbu- 
lence occurs in any case in which the Reynolds number, 
Re, is greater than w 1000, where Re = dv/v, d is the 
characteristic scale of the flow instability, v is its char- 
acteristic velocity, and v is the kinetic viscosity of the 
fluid. While there have been some suggestions that the 
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ICM may have a non-negligible viscosity (Ruszkowski et 
al. 2004; Reynolds et al. 2005), this quantity remains un- 
known because the physics of such dilute and magnetized 
plasmas is poorly constrained. In particular, even minute 
magnetic fields lead to small proton gyroradii that sup- 
press viscosity efficiently. However, it has been pointed 
out that the exponential divergence of neighboring field 
lines in a tangled magnetic field may lead to only mod- 
est suppression (e.g. Narayan & Medvedev 2001). On 
the other hand, Rcbusco et al. (2005) showed that the 
turbulent scales and velocities required to spread metals 
in the Perseus cluster closely correspond to those nec- 
essary to balance cooling, and similar results have been 
recently obtained for several other clusters (Rcbusco et 
al. 2006). This turbulence has important implications 
for particle acceleration and the mixture and transport 
of gas energy and heavy elements. 

Obscrvationally, there are several circumstantial clues 
about the nature of turbulence in the ICM, such as the 
pressure distribution in the Coma cluster (Schuecker et 
al. 2004), the lack of resonant scattering in the 6.7 iron 
Ka line in the Perseus cluster (Churazov et al. 2004) and 
the Faraday rotation map of the Hydra cluster (Vogt 
& Ensslin 2005, see also Iapichino & Niemeyer 2008). 
Turbulence has also been invoked to explain the non- 
thermal emission in clusters (e.g. Brunetti & Lazarian 
2007), and Dopplcr shifts from such bulk motions are 
likely to be directly detectable with the next generation 
of X-ray observatories, such as Constellation-X, with an 
envisaged spectral resolution of 1-2 cV (e.g. Briiggen et 
al. 2005b). 

Finally, our understanding of AGN-driven turbulence 
in cool-core clusters is complicated both by a variety 
of other possible sources of turbulence and competing 
physical effects. Clusters form through the accretion 
of smaller structures and this infall can generate tur- 
bulence (see e.g. Takizawa 2005). Episodes of active 
merging are also expected to produce turbulence, as seen 
in simulations, (e.g. Norman & Bryan 1999; Ricker & 
Sarazin 2001). Within clusters, the motion of galax- 
ies can also produce wakes that are likely to generate 
turbulent motions (Brcgman & David 1989; Kim 2007), 
and microphysical processes, often described in terms of 
conduction and viscosity, may also play important roles 
(e.g. Narayan & Medvedev 2001; Voigt & Fabian 2004; 
Ruszkowski et al. 2004; Sternberg et al. 2007). In sum- 
mary, cool-core clusters are a mystery that is carefully 
observed, poorly understood, and closely tied up with 
AGN-driven turbulence. 

It is with this in mind that we have carried out detailed 
simulations of AGN-driven turbulence in a cool-core clus- 
ter using the adaptive mesh refinement code, FLASH3. 
While the direct simulation of turbulence is extremely 
challenging, computationally expensive, and dependent 
on resolution (e.g. Glimm et al. 2001), its behavior can be 
approximated to a good degree of accuracy by adopting 
a sub-grid approach. In this case the flow is decomposed 
into mean and fluctuating components, which provides 
a systematic framework for deriving a set of turbulence 
equations. 

There are several types of models that are able to de- 
scribe such fluctuations arising from the RT instability 
as well as the Richtmyer-Meshkov instability, which oc- 
curs when a shock hits a medium of varying acoustic 



impedance. The simplest such model consists of ordi- 
nary differential equations for the mixing region (e.g. 
Alon 1995; Chen et al. 1996; Ramshaw 1998), describing 
for example, the amplitude of the bubble by balancing 
inertia, buoyancy, and drag forces. While these yield the 
right growth rates, they fail when there are multiple in- 
terfaces and are not readily extended to two and more 
dimensions. Although these problems can be addressed 
with multifluid models (Youngs 1989), such models are 
complicated, numerically expensive, and sometimes un- 
stable. 

A second class of models evolves the turbulent kinetic 
energy per unit mass and its dissipation rate. Such 
"two-equation turbulence models" developed for unsta- 
ble shear flows postulate a turbulent viscosity, a Reynolds 
stress, and a dissipation term (e.g. Llor 2003). However, 
the usual Reynolds stress terms must be modified in the 
presence of shocks, and modeling the RT and RM in- 
stabilities requires a buoyancy term that depends on the 
amplitude and the wavelength. 

Recently, DiMonte & Tipton (2006, hereafter DT06), 
described a sub-grid model that is especially suited 
to capturing the buoyancy-driven turbulent evolution 
of AGN bubbles. The model captures the self-similar 
growth of the RT and RM instabilities by augmenting 
the mean hydrodynamics equations with evolution equa- 
tions for the turbulent kinetic energy per unit mass and 
the scale length of the dominant eddies. The equations 
are based on buoyancy-drag models for RT and RM flows, 
but constructed with local parameters so that they can 
be applied to multidimensional flows with multiple ma- 
terials. The model is self-similar, conserves energy, pre- 
serves Galilean invariance, and works in the presence of 
shocks, and although it contains several unknown coef- 
ficients, these are determined by comparisons with ana- 
lytic solutions, numerical simulations, and experiments. 

Here we implement the DT06 model into FLASH3 to: 
1.) examine the impact of turbulence on the morphol- 
ogy and stability of AGN-driven bubbles as observed in 
nearby clusters (e.g. Fabian 2006); 2.) quantify tur- 
bulence, comparing it to present indirect entropy and 
metal distribution constraints (e.g. Inogamov & Sun- 
yaev 2003) and making predictions for radial profiles, 
as directly measurable by future linewidth studies. Our 
goal is to focus on better understanding the basic case 
of purely AGN-driven turbulence in an inviscid, unmag- 
netized ICMm and to this end, we follow the model of 
Roediger et al. (2007; hereafter R07), in which the ICM 
is described by a spherically-symmetric profile fit to X- 
ray observations of the Perseus cluster, and feedback is 
modeled by periodically injecting energy into the center 
of this distribution. 

The structure of this work is as follows. In §2 we give 
an overview of the FLASH3 code and our implementa- 
tion the DT06 subgrid-turbulence model within it. In 
§3 we present tests of our implementation against an- 
alytic solutions. In §4 we discuss our modeling of the 
galaxy cluster and energy input by the central AGN. In 
§5 we present the results from our simulations and dis- 
cuss their observational consequences. Our conclusions 
are summarized in §6. 
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2. NUMERICAL MODELING 
2.1. FLASH 

All simulations in this study were performed with the 
multidimensional adaptive-mesh refinement (AMR) code 
FLASH3 (Fryxell et al. 2000). FLASH3 is a modular 
block-structured code, which is parallelized using the 
Message Passing Interface (MPI) library. It advances the 
equations of inviscid hydrodynamics by solving the Rie- 
mann problem on a Cartesian grid using a directionally- 
split Piecewise-Parabolic Method (PPM) solver (Colella 
& Woodward 1984; Colella & Glaz 1984; Fryxell, Muller, 
& Arnett 1989), a higher-order version of the method de- 
veloped in Godunov (1959). This approach uses a mono- 
tonicity constraint rather than an artificial viscosity to 
control oscillations near discontinuities. However, the 
discretisation of the equations still leads to a numerical 
viscosity, as discussed in further detail below. Finally, 
FLASH3 offers the opportunity to advect mass scalars 
along with the gas density, a property that we utilize 
both in advancing the variables in our subgrid-turbulence 
model and capturing the evolution of metals. 

2.2. Turbulence Model and Numerical Implementation 

To capture the development of AGN-driven turbulence 
and its impact on the intracluster medium, we make use 
of the subgrid model developed in DT06, in which the 
Navier-Stokes equations are extended with a turbulent 
viscosity, p t , that depends on a turbulent eddy size, L, 
and a turbulent kinetic energy per unit mass, K. The 
flow is divided up into a mean component and fluctuating 
component, such that the total velocity u is given by 

u = u + u', (1) 

where u is the mass averaged mean velocity u = pu/p, 
and pu! = 0, where p is the mass density and the over- 
bar indicates an ensemble average over many realization 
of the flow. To leading order, this expansion gives the 
following evolutionary equations: 

dp | dpUj _ Q 
dt dxj 
dpui dpuiiij dP dRi.j 
dxi 



(2) 



dt dxj 
dpE dpEuj 



dxj 

_d_ (l^QE_ 

dxj \ Ne dxj 



(3) 



dt dxj dxj \ Ne dxj J dxj 



Sk, 



where t and x are time and position variables, p(x, t) is 
the average density field, u;(x, i) is the mass- averaged 
mean-flow velocity field in the i direction, P(x, i) is the 
mean pressure, and -E(x, t) is the mean internal energy 
per unit mass. Subgrid turbulence affects these mean- 
flow quantities through the explicit source term, Sk , the 
Reynolds stress tensor Rij, and the turbulent viscosity, 
which is scaled in the energy equation by Ne- In the case 
in which multiple fluids are considered, these equations 
are supplemented by a mass-fraction equation: 



dpF r dpF r Uj 

~d~r M 



dx i 



_d_ 

dxj 



tit dF r 



(4) 



"3 \ N F dXj 

where F r is the mass fraction of species r in a given zone, 
and Ne is a scale factor. 

The turbulence quantities that appear in these equa- 
tions are calculated from evolution equations for L and 



K. The eddy scale L must be evolved because the 
buoyancy-driven RT and RM instabilities depend on the 
eddy size, which is expected to grow self-similarly. Sim- 
ple equations that include diffusion, production, and 
compression are given by 



dpL dpLiij 

dt dxj 

and 

dpK dpKiij 

dt dxj 



9 ( in dL 
dx, \Nl dxj ' 



Cc-pL^r, (5) 



d 

dxj 



jM_dK_ 

Nk dxj 



dxj 



(6) 



where Nk, Nl, Cc are constants. 

In the L equation the three terms on the right hand 
represent, respectively: turbulent diffusion, growth of ed- 
dies through turbulent motions, and growth of eddies 
due to motions in the mean flow. In the K equation the 
three terms on the right hand side represent, respectively: 
turbulent diffusion, the work associated with the turbu- 
lent stress, and a source term Sk that contains both 
Rayleigh- Taylor and Richtmyer-Meshkov contributions, 
which also appears in the internal energy equation to 
conserve energy. 

The key source term for the RT and RM instabilities 
is based on the successful buoyancy-drag model, namely 



S K = pV 



CsAigi — Ci 



V 2 



(7) 



where the coefficients Cb and Co are fit to experiments. 
Physically, Cb describes turbulent entrainment, which 
reduces the density contrast, and the drag coefficient Cd 
describes the dissipation of turbulent energy when the 
average scale is characterized by L. Moreover, V = \2K 
is the average turbulent velocity, gi = —(l/p)dP/dxi 
is the gravitational acceleration, and Aj is the Atwood 
number in the i direction. In the code, we determine this 



as 



As = 



P+ ~ P- 
P++P- 



C A 



dp 



p + L\dp/dxi\ dxi 



(8) 



where Ca is a constant and p_ and p + are the densities 
on the back and front boundaries of the cell in the i 
direction. 

As the inclusion of shear-driven turbulence in the 
DM06 model is still in the process of development, for 
the Reynolds stress tensor we consider only the turbu- 
lent pressure 



(9) 



where S^j is the Kroniker delta and Cp is a constant. 
Thus, although we expect buoyancy-driven turbulence to 
dominate ICM mixing by hot bubbles, our results nev- 
ertheless represent a lower-limit that does not include 
shear-driven turbulence such as that arising from the 
Kelvin-Helmholtz instability (Helmholtz 1868; Kelvin 
1871). 

Finally, the turbulent viscosity is calculated as 



fit = C^pLV, 



(10) 



where C M is a constant. A list of the model coefficients 
and their values is given in Table [1] in which we also 
summarize the effect of each constant on the model. 

Our numerical implementation of these equations is di- 
vided into three steps that are carried out after the main 
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TABLE 1 

Summary of coefficients in the Turbulence model. 
Constants that are fit to experiment appear with error 
bars, and in those cases we take the central value for 
this study. 



Parameter 


Value 


Effect 


N L 


0.5 ±0.1 


Diffusion of L 


N K 


1.0 ±0.2 


Diffusion of K 


N F 


1.0 ±0.2 


Diffusion of Species F 


N E 


1.0 ±0.2 


Diffusion of E 


Ca 


2.0 


Turbulent Atwood Number 


C b 


0.84 ±0.11 


Buoyancy-Driven Turbulence 


C c 


1/3 


Compression of L 


C D 


1.25 ±0.4 


Drag term for K 


C P 


2/3 


Turbulent Pressure 




1.0 


Turbulent Viscosity 



hydro step in FLASH3, which advects all the variables 
above, including K and L. In the first step, we imple- 
ment the diii/dxi terms in eqs. ([3]), ([5]), and J6]) explic- 
itly. In the second step, we: (i) compute V as \/2K, (ii) 
use a leapfrog technique to add the source term to V as 
Sk/ pV along with the pV to the L equation, and then 
(hi) write V back to the K array as K = V 2 /2. Finally, 
in the third step, we calculate the turbulent viscosity and 
use this to implement the diffusive mixing terms in eqs. 
©, (O, and ^ explicitly. This final step requires us to 
impose an additional constraint on the minimum times 
step of dt < (A 2 //i t )/4 where (A) is the minimum of dx, 
dy, and dz in any given zone. This diffusive constraint 
must be satisfied for all zones in the simulation. 

3. TESTS OF OUR SUBGRID-TURBULENCE MODEL 

3.1. Rayleigh- Taylor Shock- Tube Test 

As a test of our implementation, we recreated the 
Rayleigh- Taylor test problem described in section V of 
DT06. We considered a 1 cm region that was filled with 
two 7 = 5/3 ideal fluids with constant densities: pi = 1 
g/cm 3 from x — —0.5 to x — 0, and pi = 0.9 from x = 
to x = 0.5. The region was placed in a gravitational field 
that pointed in the x direction with an acceleration of 
9.8 x 10 s cm/s 2 , or 10 6 times the Earth's gravity, and 
the temperature profile was chosen such that the overall 
distribution was in hydrostatic balance and at x = the 
temperature of the lower density fluid was T 2 — 50 K. 
Although this is a one-dimensional problem, as a check 
of our implementation our test simulations were carried 
out in 2 dimensions in a 50 "block" by 1 "block" region, 
where a FLASH3 block represents 8x8 simulation cells. 
Finally, the code was allowed to place up to 3 refine- 
ments, such that Z ro fi nc was 4, based on the density and 
pressure profile, which lead to an initial dx near the in- 
terface of 1 cm/50/8/2 4 = 3.1 x 10~ 4 cm, and a minimum 
resolution of 2.5 x 10~ 3 cm. This initial profile is shown 
in Figure Q] 

As described in DT06, the growth of turbulence at the 
interface in this Rayleigh- Taylor unstable interface has 
the following analytic solution: 

L(x,i)=L(t,0) [l-x 2 /h(t) 2 ] 1/2 , (11) 
K(x,t)=K(t,0) [l-x 2 /h(t) 2 ] , (12) 

where h(t) = ai,A(0)t 2 is a scale length for the interpen- 
etrating fluid (with a b = C A C B /[S(l + 2.C D )} = 0.06), 
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Fig. 1. — Initial set-up for our Rayleigh- Taylor test simulations. 
Top: Initial density profile in the p2 = 0.9 case (solid line), and in 
the p2 =0.8 case (dashed line). Bottom: Temperature profiles in 
both cases, with lines as in the upper panel. 

and L(t,0) = h(t)/2 and K(t,0) = (dh/dt) 2 /2 are the 
turbulent length scale and turbulent kinetic energy per 
unit mass at the interface. Following DT06, we initialize 
our simulation at a time of lOps and set L(0, 10/is) = 
h(lQps) 2 / dx in the zones at either side of the interface. 
The resulting K, L, and density profiles are shown in 
Figured 

At all times, and for all quantities of interest, our 
implementation reproduces the analytic solution with a 
high degree of accuracy. As expected, the kinetic energy 
per unit mass, turbulent viscosity, and the width of the 
turbulent layer (taken to be the scale at which K reaches 
a 1/10 of its maximum value) all increase as oc t 2 , which 
eventually leads to rapid mixing between the two materi- 
als. Furthermore, although our explicit implementation 
of turbulent diffusion requires us to impose a time step 
oc dx 2 1 pt, this works well in concert with the AMR hydro 
solver. As p t increases near the interface, density gradi- 
ents are smoothed, allowing the code to derefine these 
boundaries. Thus the diffusive time step remains greater 
than or comparable to the one required by the Courant 
condition for most of the simulation, dropping only to a 
minimum value of about 1 /6 of the Courant time step at 
the very latest times, when the entire simulation volume 
has moved to the lowest allowed refinement value. 

In Figure [3] we show the results of a similar test in 
which the density on the right hand side of the volume 
decreased to 0.8 g/cm 3 , doubling the growth rate of the 
turbulent layer. Again, we find good agreement between 
our numeric approach and the analytic results, and again 
derefinements at late times keep the the diffusive dt to 
manageably large values. We repeated this 2D simula- 
tion with the x and y coordinates interchanged, and car- 
ried out a 3D simulation in which gravity was pointed in 
the z direction. In both cases we achieved results identi- 
cal to the ones presented above. 

Finally, we carried out comparison runs in which we 
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Fig. 2. — Evolution in the p2 = 0.9 case. Top: Profiles of the 
turbulent kinetic energy per unit mass at 50 fis, 100 fis, 200 fis, 
and 300 fis. In each case the solid line gives the simulation results 
and the dashed line gives the analytic solution. Second panel: 
Profiles of the turbulent length scales from the simulation (solid) 
and the analytic solution (dashed), with times as in the upper 
panel. Third Panel: Profiles of the simulated turbulent viscosity 
per unit mass (solid) as compared to the analytic values (dashed). 
Bottom: Simulated density profiles with times as in the other 
panels. 

initialized the simulation as in the a;-direction test runs, 
increased the resolution in the y direction to 4 or 8 blocks, 
and turned off the subgrid turbulence model. In this case, 
the fluid remained stationary until perturbations from 
numerical errors grew to the point that fingers formed 
between the two materials. The delay in the onset of this 
stage was dependent on the resolution in the y direction, 
and after this onset, the region containing interpenetrat- 
ing fingers grew roughly as t 2 . 

3.2. Richtmyer-Meshkov Test 

Next, we tested the ability of our code to simulate 
Richtmyer-Meshkov (RM) amplified turbulence (Richt- 
myer 1960; Meshkov 1969), which occurs when a shock 
encounters a discontinuity in acoustic impedance such as 
a density step. In this case, the shock accelerates the in- 
terface to a velocity v m t , which amplifies the initial per- 
turbations by sending the low-density regions running 
out ahead of the higher-density regions (e.g. Mikaelian 
1989; Alon et at 1995; Dimonte 1999; Holmes et al. 
1999). As the perturbations grow and become nonlin- 
ear, the growth rate decays away in time, leading to 
a bubble amplitude that grows as h b oc (^inti)b, where 
9 b w 0.25 ± 0.05. 

To study the ability of our implementation to capture 
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Fig. 3. — Evolution in the p2 = 0.8 case. Top: Profiles of 
the turbulent kinetic energy per unit mass at times of 50 fis, 100 
/is, 150 fis, 200 fis, where again the solid lines give the simulation 
results and the dashed lines give the analytic solution. Second 
Panel: Profiles of the turbulent length scales with lines and times 
as in the upper panel. Third Panel: Profiles of the simulated and 
analytic turbulent viscosity per unit mass. Bottom: Simulated 
density profiles. 

this instability, we reproduced the shock-tube test prob- 
lem described in section VI of DT06. In this case, we 
considered a 15 cm long region, spanning from x = — 7 
to x = 8. The rightmost section of the simulation vol- 
ume, from x = to x = 8, was filled with a stationary gas 
with p — 0.667 g/cm 3 and T = 2.16 K. The center region 
of the simulation, from x = —6.9 to x — was filled with 
a second stationary gas with p — 1 and T — 1.44 K. Both 
gases has a poly tropic index of 7 = 5/3. Finally, the re- 
gion from x = —7 to x = —6.9, was filled with ap = 1.81 
g/cm 3 , T = 2.27 K, 7 = 5/3 gas flowing in from the left 
boundary at 10 4 cm/s, such that a shock with a Mach 
number of 1.57 was established in the central gas. The 
simulation volume was again two-dimensional, 60 blocks 
in the x direction by 1 block in the y direction, and again 
3 additional levels of refinement were allowed, such that 
the minimum and maximum zone size in the x direction 
were 0.0039 cm and 0.031 cm respectively. Finally, along 
the interface we initialized L = 0.01. 

For this problem the turbulent length scale and kinetic 
energy are expected to grow as 

L(t) = L(O)[tv nM /0 + l] e , (13) 
K(t) = K(0)[tv RM /6 + l] 2e - 2 , (14) 

where 9 = (2Cd + 1.5) -1 = 1/4 where wrm is the initial 
post-shock time derivative of L (DT06). In Figure 2] we 
show the K and L profiles that occur along the density 
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Fig. 4. — Top: Profiles of the turbulent kinetic energy per unit 
mass at 12 outputs spaced at 25 /is intervals, beginning at the 
moment the shock passes through the interface. For clarity, we 
alternate outputs between solid lines (25 fis, 75 fis, etc.) and dotted 
lines (50 fis, 100 fi, etc.). The dashed line shows th e analytic 
solution for the evolution of K at the interface, eq. JlD . Bottom: 
Profiles of the simulated turbulent length scale at the time steps 
given in the upp er panels. The dashed line shows the analytic 
solution, cq. 1131 . 

interface as the shock passes through, as compared to 
these analytic solutions. Once we account for the fact 
that L is initially compressed as the shock passes through 
the interface, eqs. (fl3|) and (fl4]) . with vrm = 1-6 x 10 4 
cm/s, give an excellent match to the simulation results. 
Again, we repeated this experiment interchanging the x 
and y, as well as in the z direction, obtaining identical 
results. 

4. CLUSTER MODELING 

Having carried out detailed tests of our subgrid- 
turbulence model, we moved on to modify FLASH3 to 
apply this model to the context of the AGN heated ICM. 
This involved choosing appropriate initial conditions and 
metal injection profiles, modifying the code to account 
for the transport and turbulent mixing of metals, and 
modeling AGN feedback. 

4.1. Cluster Profile and Metal Injection 

For our overall cluster profile, we adopted the model 
described in R07, which was constructed to reproduce the 
properties of the brightest X-ray cluster A426 (Perseus) 
that has been studied extensively with Chandra and 
XMM-Newton. In this case, the electron density n c 
and temperature T c profiles are based on the deprojected 
XMM-Newton data (Churazov et al. 2003; 2004) which 
are also in broad agreement with the ASCA (Allen & 
Fabian 1998), Beppo-Sax (De Grandi & Molendi 2001; 
2002) and Chandra (Schmidt et al. 2002; Sanders et al. 



2004) data. Namely: 
4.6 x 10~ 2 



4.8 x 10- 



tl + (£) 2 ] 1 - 8 + 



.^210.87 



and 



T 



7 x 



[l + (£) 3 ] 
[2-3 +(£) 3 



keV, 



(15) 



(16) 



where r is measured in kpc. Furthermore, the hydrogen 
number density was assumed to be related to the electron 
number density as nn = n e /l.2 according to standard 
cosmic abundances. The static, spherically-symmetric 
gravitational potential was set such that the cluster was 
in hydrostatic equilibrium. 

The metal injection rate in the central galaxy was as- 
sumed to be proportional to the light distribution, and 
modeled with a Hernquist profile given by 



Pmetal(?", t) 



M metal a 



1 



2ir r (r + a) 3 ' 



(17) 



with a scale radius of a = 10 kpc and a total metal in- 
jection rate of Mmotaio = 952M Q Myr~ 1 . Note that the 
injection rate could also have been time-dependent, as 
used in Rebusco et al. (2005), to account for the higher 
supernova rate in the past and the evolution of the stel- 
lar population (see also Renzini et al. 1993). However, 
as we follow the evolution of the cluster for only about 
500 Myrs, the metal injection rate does not change sig- 
nificantly over this time (Rebusco et al. 2005). 

4.2. Tracing the metals 

In order to be able to trace the metal distribution, we 
utilize a mass scalar to represent the local metal fraction 
in each cell, F = p m ctai/p- Hence, the quantity Fp gives 
the local metal density, p m ctai, which has a continuity 
equation including the metal source, given by 



dpF dpFHj 



dt 



+ 



dxj 



d 

dxj 



N F dxj 



Pmctal- 



(18) 



Furthermore, we assumed that the metal fraction is small 
at all times. Hence, we could neglect pmetai as a source 
term in the continuity equation for the gas density. 

4.3. Bubble generation 

Bubbles in the ICM are thought to be inflated by a pair 
of ambipolar jets from an AGN in the central galaxy that 
inject energy into small regions at their terminal points, 
which expand until they reach pressure equilibrium with 
the surrounding ICM (Blandford k Rees 1974). The re- 
sult is a pair of underdense, hot bubbles on opposite sides 
of the cluster center. As in R07 these were modeled using 
two methods. 

In one method, we injected a total energy of E- m j over 
an interval of length Ti n j = 5 x 10 6 years into each of two 
small spheres of radius r m j . The gas inside these spheres 
was heated and expanded similar to a Sedov explosion to 
form a pair of bubbles in a few Myrs, a time much shorter 
than the rise time of the generated bubbles. The parame- 
ters ri n j , -Einj and Tj n j were chosen such that these regions 
reached a radius of rbbi and a density contrast of approx- 
imately Pb/pamb — 0.05 as compared to the surrounding 
ICM. However, the dependence of the bubble dynamics 



Scannapieco & Briiggen 7 



TABLE 2 

Run Parameters 



"Run 

jtvun 


/ 

t rcnnc 


Maximum 


Effective 


Subgrid 


Cooling 


J3UDD1C 


JjUDDIC 


.DUD DlC 


JjUDDie 


NeLIHG 




Resolution 


Grid 


1V1U01U1 




lype 


Freq. 


Radius 


WllftC I 


CAT A TT'C 

DIN Afjo 


u 


0.6 


■6 kpc 




IN O 


NT a. 
^N O 


Evac. 


Once 


11 Kpc 


LO. Z KpC 


en A T?c 

O-DAHjD 





0.6 


■6 kpc 


i no/i 3 


1 OS 


IN O 


Evac. 


Once 


11 KpC 


lO.Z KpC 


OTN Afjo 


Q 
O 


2.6 


! kpc 


ZOO 


IN O 


IN O 


Evac. 


Once 


11 KpC 


lo.Z KpC 


OUAHjD 


Q 

o 


2.6 


■ kpc 


0^3 
zoo 


Vac 
1 Oft 


^N O 


Evac. 


Once 


11 KpC 


lO.Z KpC 


4 IN ArLo 


1 


1.3 kpc 


^1 o3 


1NO 


IV,--, 
IN O 


Evac. 


Once 


11 KpC 


lO.Z KpC 


41JA_C;D 


1 


1.3 kpc 


K~\ o3 


Vac 
1 Ob 


IN O 


Evac. 


Once 


11 KpC 


lo.Z KpC 


DIN Alio 


D 


0.33 kpc 


ori/i fi3 


"Mr* 
IN O 


IN O 


Evac. 


Once 


11 KpC 


lO.Z KpC 


^Fl A TTCJ 


D 


0.33 kpc 




Vac 
1 Oft 


IV n 
IN O 


Evac. 


Once 


1 1 KpC 


LO.Z KpC 


civr A TT'tl 
OlN rt£jo- 


c 



0.6 


■6 kpc 


1UZ4 


j.N O 


IN O 


Evac. 


Once 


I 1 L--r»a 

I I KpC 


10. V KpC 


cpi A T?Q 

DJJA-CjD- 


c 
o 


0.66 kpc 


i no/i 3 

1UZ4 


Vac 
1 Oft 


IV n 
IN O 


Evac. 


Once 


11 KpC 


lo.U KpC 


5NAES+ 


5 


0.6 


6 kpc 


1024 3 


No 


No 


Evac. 


Once 


11 kpc 


13.4 kpc 


5DAES+ 


5 


0.6 


6 kpc 


1024 3 


Yes 


No 


Evac. 


Once 


11 kpc 


13.4 kpc 


5NAES_x 


5 


0.6 


■6 kpc 


1024 3 


No 


No 


Evac. 


Once 


11 kpc 


13.0 kpc 


5DAES_x 


5 


0.6 


6 kpc 


1024 3 


Yes 


No 


Evac. 


Once 


11 kpc 


13.0 kpc 


5NASS 


5 


0.6 


6 kpc 


1024 3 


No 


No 


Sedov 


Once 


11 kpc 


13.2 kpc 


5DASS 


5 


0.6 


■6 kpc 


1024 3 


Yes 


No 


Sedov 


Once 


11 kpc 


13.2 kpc 


5NACR 


5 


0.6 


■6 kpc 


1024 3 


No 


Yes 


Evac. 


50 Myr 


16 kpc 


17.0 kpc 


5DACR 


5 


0.6 


6 kpc 


1024 3 


Yes 


Yes 


Evac. 


50 Myr 


16 kpc 


17.0 kpc 


5NCSR 


5 


0.6 


■6 kpc 


1024 3 


No 


Yes 


Sedov 


50 Myr 


16 kpc 


17.0 kpc 


5DCSR 


5 


0.6 


■6 kpc 


1024 3 


Yes 


Yes 


Sedov 


50 Myr 


16 kpc 


17.0 kpc 



on the density contrast, pb/Pamb, is weak provided that 
Pb/p&mb -C 1 (R07). In addition to the bubbles, the ex- 
plosion sets off shock waves that move through the ICM, 
which in fact are the energetically dominant component, 
as discussed in more detail below. 

In another method, we evacuated the bubble regions 
by removing gas while keeping them in pressure equi- 
librium with their surroundings. Inside two spheres of 
radius rbbi, we removed the gas at a rate p for a time 
7"cvac = 5 x 10 6 years that was small compared to the 
evolution timescale of the bubbles, but long enough to 
prevent numerical problems. The sink rate p was set to 
decrease the density inside the bubbles down to a den- 
sity contrast of /3b/Pamb compared to the surrounding 
ICM, and in order to conserve the metal mass, the metal 
fraction was also updated during the evacuation. Since 
gas was removed from a small volume in the process of 
forming the bubble, total mass is not conserved in this 
method. However, this does not affect the density pro- 
file of the cluster on the scales that we consider here. 
In fact, as illustrated in Fig. 4 of R07, the differences 
between this method and the one described above are 
small for the mean- flow variables. However, the evacu- 
ation method does not cause shocks that move through 
the ICM, making it less likely to drive further turbulence 
by the Richtmyer-Meshkov instability. Thus we used this 
model to separate out the effects of buoyancy-driven and 
shock-driven turbulence on galaxy clusters, and consider 
it first in our study. Note however that in both cases the 
bubble models are simplified and do not capture the de- 
tails of the initial heating of the ICM by AGN jets, which 
can also affect stability {e.g. Sternberg & Soker 2008b) 

4.4. Simulation Parameters 

Recent work has considered the possibility that AGN 
bubbles are stabilized by molecular viscosity and mag- 
netic fields. Reynolds et al. (2005) found that even 
a modest shear viscosity (corresponding to 1/4 of the 
Spitzer value) can quench fluid instabilities (see also 
Ruszkowski et al. 2004). Smoothed particle hydrody- 



namic simulations of buoyant bubbles in a viscous ICM 
were also performed by Sijacki & Springel (2006). Fur- 
thermore, it has been known for some time that the ICM 
hosts significant magnetic fields with typical strengths of 
5 p,G (Carilli et al. 2002). In many cases, the restor- 
ing tension generated by bending of the field lines may 
exceed the buoyancy force driving the RT instability 
(Chandrasekhar 1961). De Young (2003) derived ana- 
lytic conditions for the stabilization by ICM magnetic 
fields, suggesting that observed field strengths might sta- 
bilize the bubble interface in the linear regime, and idea 
explored in further detail by the two-dimensional MHD 
calculations in Briiggen & Kaiser (2001) and recent 3D 
simulations by Ruszkowski et al. (2007). Yet, although 
viscosity and magnetic fields can clearly affect bubble dy- 
namics and turbulence, their impact is still poorly under- 
stood, and we focus here on the basic case of turbulence 
arising in an inviscid and unmagnetized fluid. 

All our simulations are performed in a cubic three- 
dimensional region 680 kpc on a side, with all reflecting 
boundaries. For our grid, we chose a block size of 8 3 
zones and an unrefined root grid with 8 3 blocks, for a 
native resolution of 10.6 kpc. The refinement criteria 
are the standard density and pressure criteria, and we 
allow for 4 levels of refinement beyond the base grid, 
corresponding to a minimum cell size of 0.66 kpc, and 
an effective grid of 1024 3 zones. The parameters of all 
the simulations carried out in this study are summarized 
in Table [H We name each run by concatenating the 
maximum level of refinement in FLASH3, the presence of 
a subgrid-turbulence model (N for none or D for DT06), 
whether the run uses the cooling model described below 
(A for adiabatic or C for cooling), the choice of bubble 
type (S for Sedov- Type and E for evacuated), and the 
whether the runs contain a single pair of bubbles (S) or 
periodic (P) AGN feedback. 
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Fig. 5. — Snapshots of mean-flow quantities in evacuated single-bubble runs at t = 50 Myrs, 100 Myrs, 150 Myrs, and 200 Myrs 
(arranged from top to bottom in each column) . All panels show values at a z = slice through our simulations and cover the region from 

r. ~ V9 _ 

g cm 



-100 to 100 kpc and y = —100 to 100 kpc. First Column: Contours of logp spanning the range from p = 10 27 g cm 3 to 10 



cm 3 from the 5NAES pure-hydro run. Second Column: Contours of log p from the 5DAES run with subgrid turbulence, spanning the 



same range of densities as in the first column. Third Column: Contours of log T from T - 
Column: logT contours from the 5DAES run with the same scale as in column 3. 



10 7 5 K to 10 9 K from the 5NAES run. Fourth 



5. RESULTS 

5.1. Evacuated Bubbles 

As our first case study, we carried out two simulations: 
an adiabatic pure-hydro simulation (5NAES) and an adi- 
abatic simulation with subgrid turbulence (5DAES). At 
the start of each simulation a single pair of nibble = 11 
kpc bubbles was centered along the x — y axis at an offset 
of Rq = ±13.2 kpc, and within the bubbles the gas was 
evacuated to Pb/pamb — 0.05 and heated to 20 times the 



temperature of the surrounding medium. Snapshots of 
density and temperature from these simulations appear 
in Figure [3 spanning times from 50 to 200 Myrs. 

The differences between the two runs are dramatic. 
In the pure-hydro run the bubbles fragment after rising 
a single pressure scale height. The dominant unstable 
modes that are set by the resolution of the adaptive grid 
quickly shred the evacuated regions, drastically reduc- 
ing them in size. If we were able to repeat this sim- 
ulation with arbitrarily high resolution, we would find 
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Fig. 6. — Snapshots of properties of subgrid turbulence in the 5DAES run at t = 50 Myrs, 100 Myrs, 150 Myrs, 200 Myrs (from top to 
bottom in each column). As in Figure [5] all panels show a central 200 kpc X 200 kpc 2 = slice. First Column: Logarithmic contours of 
the turbulent velocity labeled in cm s _1 from V = \J1K = 0.3 km s _1 to V = 300 km s — Second Column: Logarithmic contours of the 
local turbulent Mach number V/c a , from 10 -4 to 10 _1 . Third Column: Logarithmic contours of L, labeled in cm, from 0.1 kpc to 10 kpc. 
Fourth Row: Logarithmic turbulent viscosity per unit density jLit/p, labeled in cm 2 s — 1 from 10 — 2 km/s kpc to 10 3 km/s kpc. 



that the bubbles have developed fingers upon fingers of 
turbulence, which penetrate deep into the surrounding 
medium. However, this dispersal is halted by the finite 
resolution in the pure-hydro run, which by 100 Myrs has 
already significantly underestimated the mixing between 
the evacuated region and its surroundings. 

On the other hand, the subgrid-turbulence run cap- 
tures this mixing, modeling the growth of all the modes 
that our computational grid cannot resolve. The super- 
position of these modes then smears out the interface be- 
tween the bubble and the ambient medium, transforming 



the bubbles into clouds of mixed material, which stay in- 
tact and expand as they rise in the stratified ICM. Note 
that in both runs the bubbles are trailed by streams of 
colder gas, which are lifted from the center and cool adi- 
abatically, consistent with observations (e.g. Simionescu 
et d. 2007). 

To quantify the properties of the turbulent clouds fur- 
ther, in Figure [6] we plot slices from our 5DAES sim- 
ulation, focusing on quantities evolved by the subgrid 
model. Turbulent velocities quickly approach a maxi- 
mum of w 100 km s _1 , roughly 5% of the sound speed 



10 



EI' 






♦ 



Fig. 7. — Snapshots of metal and emission distribution in the central 200 X 200 kpc slice at z = 0, at t = 50 Myrs, 100 Myrs, 150 Myrs, 
and 200 Myrs (from top to bottom in each column) with metallicity increasing inwards. First Column: Contours of logp meta i from the 
5NAES run, at a fixed but arbitrary scale, with metallicity increasing inwards. Second Column: Contours of logp meta i from the turbulence 
run, with a scale matching the first column. Third Column: Unsharp X-ray image from the 5NAES run, projected in the z-direction. 
Fourth Column: Unsharp X-ray image from the 5DAES run, projected in the z-direction. Fifth Column: Unsharp X-ray image from the 
5NAES run, projected in the x-direction. Sixth Column: Unsharp X-ray image from the 5DAES run, projected in the z-direction. 



within the bubbles and 10% of the sound speed of the 
surrounding ICM. As the hot gas rises at w 500 km/s, 
the turbulent motions are fast enough to mix this ma- 
terial with the surrounding medium, but much too slow 
to halt its overall motion. In addition to the velocity of 
turbulent motions, the efficiency of this mixing process 
is dependent on the scale of the turbulent eddies, which 
is plotted in the third column of Figure 6. Here we see 
that L quickly rises to w 10 kpc, but does not exceed the 
overall size of the clouds, keeping material well-mixed 
within these w 30 kpc regions. Together L and V act to 
produce a typical turbulent viscosity of w 300 km/s kpc, 
which diffuses material between the clouds and their sur- 
roundings. Thus rather than being a destructive process, 
turbulence acts as an effective mixing mechanism, which 
alters the rising structure, but operates at scales that are 
too small and with velocities that are too slow to disrupt 
it completely. 

To study this mixing process further, we compare the 
distribution of ICM metals from our 5DAES (subgrid- 
turbulence) and 5NAES (pure-hydro) simulations in the 
left two columns of Figure [7l In both runs, highly- 
enriched material from the center of the cluster is carried 
outward by the rising bubbles, but this transport is much 
more effective in the 5DAES run. In the pure-hydro run 
metal transport is halted at sa 50 kpc by bubble dis- 



ruption, which occurs when RT instabilities shred the 
evacuated region into resolution-limited cavities. In the 
subgrid-turbulence run, on the other hand, small-scale 
fluctuations act to move metals into the cloud while at 
the same time keeping the structure coherent, such that 
it continues to rise out to substantially larger radii. 

We quantify this in the upper row in Figure [5J which 
shows the overall density of metals as a function of radius 
in each of our runs. By 150 Myrs, the bubbles in the 
5NAES run have dissipated, leaving behind substantial 
enhancement of metals within 50 kpc, but having little 
or no impact outside of this radius. On the other hand, 
clouds in the subgrid-turbulence run continue to rise even 
at 200 Myrs, leading to enhancements of by over a factor 
of 5 at 70 kpc. 

In the lower panel of Figure [5J we plot Sturb) the to- 
tal turbulent kinetic energy within a given radius in our 
5DAES simulation. This can be directly compared with 
the total energy available from the buoyant rise of the 
bubbles, which can be simply estimated as (e.g. Nulsen 
et al. 2006) 



Eh 



bouyancy 



Ho 



3p 

cvac ^cvac 



(19) 



where R is the distance of the bubble from the center of 
the cluster which is initially at Rq. This gives a value of 
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Fig. 8. — Metal density and total turbulent kinetic energy at times of 50 Myrs, 100 Myrs, and 150 Myrs, 200 Myrs arranged in columns 
from left to right. Top: Metal density from the evacuated-bubble subgrid-turbulence run (5DAES, solid lines), and from the evacuated 
bubble pure-hydro run (5NAES, dashed lines). For comparison the dotted line shows the Hernquist profiles in which metal are ejected into 
our simulations. The logarithmic scale in this row is fixed, but arbitrary. Bottom: Total turbulent kinetic energy contained with in a given 
radius, as a function of that radius (solid lines). The shaded regions are bounded by 1% to 2% of i^bouyancy as given by eq. H19I I. 



ss 10 59 2 ergs per pair of rbbi = H kpc bubbles. Compar- 
ing this value to the turbulent kinetic energy in our sim- 
ulations, we find that -Eturb is only rs 1% of this overall 
energy, thus we do not expect the energy of the turbulent 
motion themselves to play a major role in the heating of 
the cool-core region in the ICM. Rather, the key impact 
of turbulence is to increase the efficiency with which the 
thermal energy of the rising cloud is mixed into its sur- 
roundings {e.g. Soker 2004; Sternberg et al. 2007; Stern- 
berg & Soker 2008b). In this sense it behaves much more 
like heat conduction (e.g. Narayan & Medvedev 2001; 
Voigt & Fabian 2004; Ruszkowski et al. 2004; Sternberg 
et al. 2007), than an energy source. Note however, that 
our simple turbulence model does not account for shear- 
driven effects, and thus may somewhat underestimate 
^turb, though not by orders of magnitude. 

Finally, the physical differences between the two runs 
lead to different observed morphologies. In the right pan- 
els of Figure[7]we present approximate X-ray images from 
each of the two runs, calculated by projecting the emis- 
sivity, computed as 



A(T)n 2 e , 



(20) 



where we estimate the cooling function, A(T), which de- 
scribes radiative losses from the optically-thin plasma, as 
in Sarazin (1986; see also Raymond et al. 1976; Peres et 
al. 1982). Furthermore, to draw out structure we create 
unsharp mask images by dividing the map by a 30 kpc 
smoothed version of itself. 
These images are comparable to those in the study by 



Reynolds et al. (2005), which was carried out with similar 
resolution, but in a hydrostatic atmosphere falling off as 
p(r) oc [1 + (r/ro) 2 ] -0 ' 75 , which is significantly less steep 
than the one we have prescribed in the central cooling 
flow region of our simulation. Furthermore our times 
can easily be related to their dimensionless units as t = 
t x 57kpc/c s w t 60 Myrs. Our 5NAES run thus confirms 
the Reynolds et al. (2005) result that pure-hydro inviscid 
bubbles fall apart within a dimensionless time of t = 3, or 
« 200 Myrs, although our steeper radial slope combined 
with our processing of the image is able to draw out 
more structure than visible in their plots. It is clear 
that these collections of fragments look nothing like the 
smooth and detached cavities of the type observed in 
Perseus (e.g. Fabian et al. 2006) and other clusters (e.g. 
McNamara et al. 2001; Choi et al. 2004; Reynolds et 
al. 2008). On the other hand, in the subgrid-turbulence 
run, the hot clouds lead naturally to coherent holes in 
the X-ray distribution not unlike those observed in cool- 
core clusters. If anything the bubbles are too large and 
radially extended as compared to observations, an issue 
we take up below, after first addressing the reliability of 
our numerical approach. 

5.2. Dependence on Resolution and Initial Parameters 

An important issue in any simulation is the degree to 
which its results are modified when one changes the min- 
imum zone size. This is of particular concern here, be- 
cause the DT06 model must evolve all subgrid modes 
contributing to the RT instability, while still achieving 
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FlG. 9. — Snapshots of logp in single-bubble runs of varying resolution. Plots show a z = slice of the region from x = —3 kpc to 100 
kpc and and y = —3 kpc to 100 kpc from the pure-hydro run. From left to right, each column shows results from simulations with 3 levels 
of refinement (2.6 kpc resolution), 4 levels of refinement (1.3 kpc resolution), the fiducial 5 levels of refinement (0.66 kpc resolution), and 
6 levels of refinement (0.33 kpc resolution). Top Row: Outputs at 50 Myrs from runs without the subgrid-turbulence model. Second Row: 
Outputs at 50 Myrs from runs including subgrid turbulence. Third Row: Outputs at 100 Myrs from runs without subgrid turbulence. 
Fourth Row: Subgrid-turbulence run outputs at 100 Myrs. Fifth Row: Outputs at 150 Myrs from runs without subgrid turbulence. Sixth 
Row: Subgrid-turbulence run outputs at 150 Myrs. 
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Fig. 10. — Snapshots of the turbulent viscosity per unit density, log /it / p, in single-bubble runs of varying resolution. As in Figurc|9] plots 
show a z = slice of the region from x = — 3 kpc to 100 kpc in the x and y directions, and from left to right, each column shows results 
from simulations with i rc fl nc = 3 (2.6 kpc resolution), £ r efmc = 4 (1.3 kpc resolution), the fiducial / r cfinc = 5 case (0.66 kpc resolution), and 
'refine = 6 (0.33 kpc resolution). From top to bottom, each row shows outputs at 50 Myrs, 100 Myrs, and 150 Myrs respectively. In all 
panels pt/p is labled in units of cm 2 s — 1 and goes from from 10 -2 km/s kpc to 10 3 km/s kpc. 



results that are resolution-independent. With this in 
mind, we carried out a set of convergence tests, repeating 
our simulations with and without the subgrid model for 
cases in which we varied the maximum level of refinement 
from Irefine = 3 (2.6 kpc resolution, 256 3 effective grid) 
up to Z rc fi nc = 6 (0.33 kpc, 2048 3 effective grid). Snap- 
shots from these simulations are presented in Figure O 
in which we have zoomed into the area from x — — 3 kpc 
to 100 kpc, and y = —3 kpc to 100 kpc to emphasize the 
regions that are most affected by turbulence. 

Focusing first on the pure-hydro runs at 50 Myrs, one 
sees a clear increase in small-scale structure as the runs 
progress to higher resolution. The wavelength of the 
fastest growing mode from a linear analysis of the RT 
instability is given by (Chandrasekhar 1981) 



(21) 



where A is the Atwood number, g is the gravitational 
acceleration, and v is the viscosity of the fluid, which, in 
the pure-hydro run, is given by the effective numerical 
viscosity. In the PPM method used by FLASH, the dis- 
sipative processes that act as an effective v are extremely 



nonlinear. Not only are the error terms in this high-order 
scheme extremely complicated, but they change with the 
nature of the flow, as PPM uses several switches that de- 
tect discontinuities (Woodward & Colella 1984; Porter 
& Woodward 1994). Thus the effective v value in our 
simulations is dependent on the structures we are trying 
to resolve, and is best expressed as an effective Reynolds 
number Re — dv/u, where d and v are the size and the 
velocity of the structure of interest. 

For our fiducial simulations (with / rc finc = 5) the ef- 
fective number of zones is 1024 3 , which corresponds to a 
spatial resolution of 0.66 kpc. In the case of large sub- 
sonic eddies, not unlike our bubbles, detailed studies have 
shown the effective Reynolds number achievable in PPM 
simulations is proportional to the number of grid points 
across the fluctuation of interest to the power of N = 3, 
where N is the order of the numerical scheme (Porter & 
Woodward 1994; Balbus et al. 1996; Sytine et al. 2000). 
In their study of viscous bubbles, Reynolds et al. (2005) 
carried out a series of runs in which numerical viscos- 
ity was explicitly added to the system at different levels, 
and they determined that FLASH3 behaved roughly as 
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Fig. 11. — Snapshots of logp in single-bubble runs with slightly different initial conditions. As in Figure [9] plots show a z = slice of 
the region from x = —3 kpc to 100 kpc and and y = —3 kpc to 100 kpc. From left to right, each row shows results from simulations in 
which the bubbles have been initially moved radially by —1%, 0%, 1%, and respectively. Finally, the rightmost row show a run in which 
the bubbles are placed along the y axis. Top Row: Outputs at 50 Myrs from runs without the subgrid-turbulcnce model. Second Row: 
Outputs at 50 Myrs from runs including subgrid turbulence. Third Row: Outputs at 100 Myrs from runs without turbulence. Fourth Row: 
Turbulence run outputs at 100 Myrs. Fifth Row: Outputs at 150 Myrs from runs without turbulence. Sixth Row: Turbulence run outputs 
at 150 Myrs. 
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if Re w 2000 — 5000. Using this number as a rough guide 
in our fiducial case, which again has resolution similar 
to that in Reynolds et al. (2005), this implies that the 
numerical viscosity is of the order of 



v ~ dv/Ke ~ 3kms 1 kpc, 



(22) 



where we assumed d ~ 25 kpc and v « 500 km/s. Taking 
this i/ value, along with j4 p» 1 and 5 » 10 -8 cm s -1 , this 
yields a fastest growing mode wavelength of A roax » 2 
kpc, which roughly corresponds to the scale of the small 
features seen in the fiducial pure-hydro run at 50 Myrs. 
From eq. (|2~Tj) . and assuming that viscosity goes as the 



oc I 



refill 



number of grid points to the power of 3, A n 
Thus decreasing the resolution to i re finc — 4, which corre- 
sponds to a spatial resolution of 1.3 kpc, moves A max up 
to w 8 kpc, greatly reducing the perturbations seen at 50 
Myrs. Similarly, choosing Z rc fme = 3 moves A max up to 
the scale of the bubbles, and no perturbations are visible. 
On the other hand, in the define = 6 run, A max is reduced 
almost to the resolution limit, leading to perturbations 
on all simulated scales. At later times, the differences 
between runs persist, such that the bubbles are shredded 
into subclumps with sizes that are strongly resolution- 
dependent. Note however, that the fastest growing mode 
as given by eq. (|2ip . need not be precisely the mode that 
dominates the late-time nonlinear evolution of the sub- 
clumps, as the growth of perturbations at A max may sat- 
urate. 

Moving to the subgrid model runs, we find no such de- 
pendence on resolution. In this case, the code carries two 
additional quantities, namely the turbulent length scale, 
L, and kinetic energy, K. It evolves L and K to repro- 
duce the analytic growth of the RT and RM perturba- 
tions in the limit in which the molecular viscosity is small 
and does not play a role on the scales of interest. The 
growth of small-scale perturbations is then completely 
captured by the evolution of L and K 1 which are used 
to construct a turbulent viscosity that is imposed on the 
mean-field variables explicitly. The role of the numerical 
viscosity in eq. (|2ip is then played by [it I Pi which grows 
to a value of 300 km/s kpc, larger than the numerical vis- 
cosity in all of the runs. Thus rather than breaking up 
into a collection of resolution-dependent subgroups, the 
mean-field density distributions remains coherent, while 
at the same time the simulation is able to maintain the 
information necessary to calculate the mixing of the tur- 
bulent material with the surrounding medium. 

Furthermore, the evolution of the turbulent viscosity 
is determined by the overall properties of the flow, and 
is largely independent of resolution. A test of this con- 
vergence is shown in Figure I10[ which gives snapshots 
of log fit I p in each of our single-bubble runs with vary- 
ing resolution. Here we see that pt/p, whose evolution 
is determined by the gravitational acceleration and the 
Atwood number, evolves to w 300 km/s kpc regardless of 
our choice of /refine- Thus the overall structure of evolving 
bubbles is quite similar across runs. 

The subgrid-turbulence model also greatly reduces the 
sensitivity of our results on the detailed choice of initial 
conditions. Figure [11] shows the results of pure-hydro 
and subgrid simulations in which we have taken the ini- 
tial offset of the bubbles to be i?o = 13.0, 13.2, and 13.4 
kpc, as well as a pair of runs in which we have taken 
Ro = 13.0 kpc, but oriented along the y-axis. In the 



pure-hydro case, our results show a strong sensitivity to 
these ±2% changes in initial position. This is because 
the code is attempting to sample an underlying turbu- 
lent density field with resolution-dependent subclumps, 
whose locations vary significantly depending on the non- 
linear, position- and velocity-dependent effective numer- 
ical viscosity. In the subgrid runs however, the presence 
of the turbulent viscosity leads to a smoothing over the 
unresolved turbulent density distribution, thus captur- 
ing the full range of positions over which fully-resolved 
turbulent material would be spread. 

5.3. Sedov Bubbles 

Next we turned our attention to the more 
computationally-demanding case in which the bub- 
bles are modeled as initially overpressured spheres. In 
this case we raised the internal temperature such that 
each bubble expanded to a density of Pb/p&mb = 0.05 
before reaching pressure equilibrium with a bubble size 
equal to the one in the evacuated model. For a 7 = 5/3 
gas, the energy input during this expansion can be 
simply calculated as 



E, 



expand 



dVp 



3pfinalVc^ 



Vin 



2/3 



1 



(23) 

where p evac and V eva _ c are the pressure and volume at the 
evacuated stage at the end of the expansion and Vinit is 
the initial volume. Although the bubble expansion oc- 
curs at approximately the sound speed of the internal 
material, this velocity is well above to the sound speed 
of the exterior gas. This energy is then deposited into 
shocks which act both to directly raise the temperature 
in the surrounding medium, as well as increase the turbu- 
lent kinetic energy through the Richtmyer-Meshkov in- 
stability as the shocks meet acoustic discontinuities. This 
energy input can be contrasted to the energy input from 
the buoyant motions of the bubbles studied in §5.1, as 
given by eq. (fl"9|). For our choice of pb/pamh = 0.05, 

(Kvac/Vinit) 2 ^ 3 — 1 = 6.35. Thus £cxpand/-Ebuoyaiicy ~ 6, 

such that overpressured bubbles should have a signifi- 
cantly greater impact on the cluster than the evacuated 
model. 

The impact of this energy is seen in the outer oval- 
shaped shock visible in the 50 Myr plots of density, 
temperature, and simulated X-rays shown in Figure 1121 
Trailing behind this discontinuity is a second density and 
temperature jump, which is due to initially inward- facing 
shocks that reflect off each other and move into the sur- 
rounding medium. The motion of these reflected shocks 
through the rising bubbles has a noticeable effect on the 
bubbles' morphology, substantially compressing them in 
the radial direction. At later times in the pure-hydro run, 
the result is a flattened distribution of distinct cavities, 
which separate and shred as the bubbles move outwards, 
much as they did in the 5NAES simulation. 

In the subgrid-turbulence case, however, the result is 
a flattened but coherent cloud which moves outwards, 
mixing with the surrounding medium while remaining in- 
tact. This distinctive shape immediately brings to mind 
familiar images of the "mushroom cloud" structure that 
results from a large airburst in the atmosphere of Earth. 
Indeed, the language used to describe such clouds, such 
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Fig. 12. — Snapshots of mean-flow quantities in Sedov single-bubble runs at 50 Myrs, 100 Myrs, 150 Myrs, 200 Myrs (arranged from 
top to bottom in each column). All panels show values at a z = 0, slice through our simulations and cover the region from x = —100 to 
100 kpc and y = —100 to 100 kpc, with scales chosen as in previous plots. First Column: Contours of logp from the run without subgrid 



turbulence, spanning the range from p 



io- 



g cm 3 to 10 25 g cm 3 . Second Column: Contours of logp from a run with subgrid 



turbulence, and spanning the same range of densities as in the upper panels. Third Column: Contours of logT from T = 10 7 ' 5 K to 10 9 K, 
from the run without subgrid turbulence. Fourth Column: log T contours from the turbulence run, again with the same dynamical range. 
Fifth Column: Unsharp X-ray image from the 5NASS run, projected in the z-direction. Sixth Column: Unsharp X-ray image from the 
5DASS run, projected in the z-direction. 



as a "toroidal fireball" perched atop a somewhat cooler 
updrafting "stem," (e.g. Glasstone & Dolan 1977), trans- 
lates naturally to describe the structures we see in our 
simulations. In fact, the physics of these two phenomena 
are deeply connected. In both cases an unstable bub- 
ble is rising upward in a hydrostatic medium. In both 
cases this bubble is neither stabilized by magnetic fields 
nor through viscosity. And in both cases the hot gas 
is well-mixed with the surrounding medium as it moves 
through it, yet it manages to retain its shape. Further- 
more, as shown in the rightmost panels of Figure 1121 
these radially-flattened clouds lead naturally to X-ray 
"holes" of precisely the type observed in cool-core clus- 
ters, without requiring additional physics. 

In Figure [13] we show the evolution of the turbulence 
variables in the 5DASS run. In general, the turbulent 
velocities and length scales are systematically enhanced 
with respect to the evacuated-bubble run, due to the ef- 
fect of the RM instability. In this case, turbulent veloc- 
ities can approach 300 km/s, 50% greater than 5DAES 
run, but still well below the ~ 500 km/s radial velocity 
of the cloud. Similarly, the turbulent length scales grow 
faster at early times and remain somewhat higher at later 



times, as compared to the evacuated case, but they never 
exceed the overall scale of the rising cloud. Together V 
and L act to generate a typical turbulent viscosity of 
10 3 km/s kpc, a substantial value that is likely to signif- 
icantly influence the overall evolution of the intraclustcr 
medium. 

In Figure [14] we quantify the radial distribution of 
metals and turbulent energy in our Sedov-bubble sim- 
ulations. In this case the differences in /o me tai between 
the pure-hydro and subgrid-turbulence runs are much 
less dramatic than they were in the evacuated-bubble 
models. The overall shock heating of the central region 
drives metals outwards to 80 kpc in 200 Myrs, re- 
gardless of whether they are transported in well-mixed 
clouds or fragmented cavities. A comparison of -Eturb to 
the energy from the buoyant bubbles shows that while 
the Richtmyer-Meshkov instability increases turbulence 
somewhat from the evacuated case, the magnitude of 
this change is much smaller than the factor of s=s 7 in- 
crease in available energy. Most of the excess energy 
instead goes directly into post-shock heating or estab- 
lishing large-scale motions. Thus even more than in the 
evacuated case, turbulence plays a role primarily in de- 
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Fig. 13. — Snapshots of properties of subgrid turbulence in the Sedov single-bubble (5DASS) run at t = 50 Myrs, 100 Myrs, 150 Myrs, 
200 Myrs (from top to bottom in each column). All panels show a central z = 0, 200 kpc X 200 kpc slice. First Column: Logarithmic 
contours of the turbulent velocity from 0.1 km s — 1 to 100 km s — 1 . Second Column: Logarithmic contours of the local turbulent Mach 
number from 10 -4 to 10 — 1 Third Column: Logarithmic contours of L from 0.1 kpc to 10 kpc. Fourth Row: Logarithmic turbulent viscosity 
per unit density from 10 — 2 km/s kpc to 10 3 km/s kpc. 



tcrmining mixing and structure, rather than in providing 
a source of energy to balance cooling. 

5.4. Periodic Bubbles 

Finally, we carried out a set of simulations to study the 
impact of our subgrid-turbulence models on the evolution 
of cool-core clusters on longer time scales. This involved 
two major changes from our single-bubble simulations. 

First, we implemented a model in which AGN-driven 
bubbles occur periodically, again following an approach 
developed in R07. As in our previous runs, we generated 



bubbles in pairs, but now at regular intervals of Tbbi = 
50 Myrs. Furthermore, we rotated the position of the 
center of the bubble pairs by 90 degrees around the y- 
axis. Thus, the first pair of bubbles was centered at 
z = and x = y = ±Rq/V2, the second pair of bubbles 
at x — and y — z — ±Rq/^/2, the third at z = at 
x = — y = ±Rq/\/2, etc., cycling back to the original 
position every 200 Myrs. 

The second major change to our simulations was the 
addition of cooling, which was calculated in the optically- 
thin limit from the same emissivity as used to construct 
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Fig. 14. — Metal density and turbulent kinetic energy (bottom row) from the Sedov single-bubbles runs at times 50 Myrs, 100 Myrs, and 
150 Myrs, 200 Myrs arranged in columns from left to right. Top: As in Figure [8] the solid lines are drawn from the subgrid-turbulence 
run (5DASS), the dashed lines are from the pure- hydro run (5NASS), and the dotted line is the profile in which metals are injected. 
Bottom: Total turbulent kinetic energy within a given radius in our simulations (solid lines) contrasted with 1-2% of Ebouyancy an d 1-2% 
of -Eexpand + ^bouyancy, represented by the lower and upper shaded regions, respectively. 
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Fig. 15. — Snapshots of mean-flow quantities in the periodic evacuated-bubble runs at 150 Myrs, 300 Myrs, and 450 Myrs (arranged from 
top to bottom in each column). All panels show values at a z = slice through our simulations and cover the region from x = —100 to 
100 kpc and y = —100 to 100 kpc, with scales chosen as in previous plots. First Column: Contours of logp from the run without subgrid 



turbulence, spanning the range from p 



10" 



g cm 3 to 10 25 g cm 3 . Second Column: Contours of logp from a run with subgrid 



turbulence, and spanning the same range of densities as in the upper panels. Third Column: Unsharp mask X-ray image from the 5NCER 
run, projected in the 2-direction. Fourth Column: Unsharp mask X-ray image from the 5DCER run, projected in the z-direction. 
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our X-ray images, namely e = A{T)n 2 e: where the cooling 
function A(T) was taken from Sarazin (1986). However, 
in order to avoid drastic cooling within a single timestep 
we did not cool cells in which the density was above 5 x 
10~ 25 g cm~ 3 or in which the temperature was below 
5 x 10 6 K. This placed an absolute entropy floor in our 
simulations at 5 keV cm 2 . 

With this cooling in place, we found that our simula- 
tion was initially radiating at w 4 x 10 44 ergs per second 
or w 1.2 x 10 60 ergs per 100 Myrs. Thus we chose to 
scale up our bubbles to a radius of rbbi = 16 kpc such 
that -Ebuoyancy ~ 5 x 10 59 ergs, and the energy input from 
buoyant heating roughly balanced cooling. As before, the 
parameters for our runs are summarized in Tabled 

Figure [15] shows snapshots of these simulations in the 
evacuated-bubble case at times of 150, 300, and 450 
Myrs. Like the single-bubbles runs, the pure-hydro simu- 
lation (5NCER) shows significantly more fragmentation 
than the subgrid turbulence run (5DCER), particularly 
towards the outskirts of the cluster. However, unlike the 
single-bubble runs, the properties of the two simulations 
continue to diverge with time at all radii. This is because 
turbulence driven by previous feedback events remains 
behind, enhancing the mixing of subsequent bubbles. 

Thus, at late times, the X-ray maps from the subgrid- 
turbulence run contain two prominent holes, while the 
distribution in the pure-hydro run is much more frag- 
mented. Note also that by 450 Myrs, errors introduced 
by the split-hydro solver employed by FLASH3 have 
grown to the point that substantial assymetries are seen 
in the pure-hydro runs. As we saw in Figure [TT] this 
dependence on small fluctuations is reduced in the sub- 
grid turbulence run, as the turbulence viscosity smoothes 
over what are essentially different realizations of the same 
turbulent density field, a field that has been unnatu- 
rally quantized into resolution-dependent subclumps in 
the pure-hydro runs. 

Figure [TH] shows the evolution of the radial profiles of 
mean-field quantities in our simulations. Focusing first 
on temperature, clear differences are apparent between 
the two runs, and these differences grow stronger with 
time. Within the central 20 kpc the temperature in the 
subgrid-turbulence run is roughly a factor of two greater 
than in the pure-hydro run, and a similar difference is ap- 
parent in the radial entropy profile. As suggested by the 
single-bubble runs, and as quantified for the 5DCER run 
below, these differences are not due to additional energy 
input from dissipated turbulent kinetic energy. Rather, 
the primary role of turbulence is to improve the mix- 
ing of the thermal energy from the interior of the heated 
cavities into the surrounding cold gas {e.g. Soker 2004; 
Sternberg et al. 2007). This means that while in the 
pure-hydro run, the bubbles deposit most of their inter- 
nal energy at the resolution-dependent radius at which 
they get disrupted, as indicated by the spike in tem- 
perature at R « 30 kpc. Hence, the turbulent motions 
captured by the subgrid model cause much more gradual 
heating, which leads to shallow temperature and entropy 
profiles at all radii < 30 kpc. 

At the same time, the well-mixed clouds in the subgrid- 
turbulence run remain coherent to larger radii, increasing 
the sphere of influence over which AGN-driven heating 
is effective. This can be clearly seen in a comparison 
of the metal profiles between the two runs, which indi- 
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Fig. 16. — Radial profile of mean-flow quantities in evacuated 
periodic bubble runs at times of 150 Myrs, 300 Myrs, and 450 Myrs, 
arranged in columns from left to right. From top to bottom, the 
rows show the average temperature, entropy, and metal density. In 
all panels the solid lines are taken from the subgrid-turbulence run 
(5DCER) and the dashed lines are taken from the pure-hydro run 
(5NCER). The dotted lines in the bottom panels show the metal 
profile that would be found if there were no resolved or subgrid 
velocities in the simulations. 

cates enhancement of almost an order of magnitude in 
the subgrid-turbulence runs relative to the pure-hydro 
runs, an enhancement which is particularly strong at 
large radii, but significant even in the central regions of 
the cluster. 

Turning to the Sedov-bubble case, similar trends are 
apparent. Here, an extremely large number of AMR 
zones were necessary at late times (« 10 5 blocks of 8 3 
zones each), forcing us to end our simulations at 300 
Myrs. However, the density slices and X-ray images pre- 
sented in Figure H7I over this time period paint a similar 
picture to the evacuated-bubble runs. Again the build- 
up of turbulent motions over time causes the properties 
of the two simulations to diverge. Again this increase in 
mixing is especially notable near the core, where a large 
number of small fragments in the pure-hydro case are 
replaced by an overall smooth and quadrapole-like con- 
figuration in the subgrid-turbulence run. And again the 
subgrid model serves to damp out much of the chaotic 
amplification of small assymetries. Furthermore, the X- 
ray images of the subgrid turbulence runs show larger 
and smoother cavities, although interestingly, the pure- 
hydro runs contain significant associations of shredded 
bubble material that appear as coherent X-ray depres- 
sions, albeit somewhat more ragged depressions than are 
seen in the observations. 

As in the evacuated-bubble case, the radial profiles 
plotted in Figure [TH] show that subgrid turbulence leads 
to more efficient heating in the central R < 30 kpc 
core. On the other hand, at larger radii heating is 
dominated by shocks and is similar between the two 
runs. Note that as S expane j is substantially greater than 
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Fig. 17. — Snapshots of mean-flow quantities in the periodic Scdo v b ubble runs at 100 Myrs, 200 Myrs, and 300 Myrs (arranged from 
top to bottom in each column). Panels and symbols are as in Figure [15] 



-EWoyancy ~ ©cool, these simulations show an overall in- 
crease in the central temperature and entropy of the clus- 
ter at late times. In reality, entropy is unlikely to rise 
beyond the point at which the cooling near the core be- 
comes inefficient, as this would halt gas accretion onto 
the central black hole and the further creation of bub- 
bles. Accounting for this change in gas accretion rate 
would lead to a self-regulating configuration (e.g. Voit 
& Bryan 2001; Birzan et al. 2004; McNamara & Nulsen 
2007; Voit et al. 2008) whose modeling falls beyond the 
scope of our study here. 

In the lower panels of this figure, we plot the metal 
profiles between the two runs. Although metal mixing is 
somewhat more efficient at larger radial in the subgrid- 
turbulence run, the metal profiles are much more sim- 
ilar between the two runs than they were in the the 
evacuated-bubble case. 

Finally, in Figure [T5] we compare the subgrid turbu- 
lence quantities between the evacuated repeated bub- 
ble run (5DCER) and the Sedov repeated bubble run 
(5NCER). As in the single bubble runs turbulent veloc- 
ities are only slightly increased by the RM instability, 
even though the total energy input in this case is sub- 
stantially higher. Thus Vturb ~ 100 — 200 km/s in both 
runs, roughly 15 — 30% of the sound speed of the cluster. 
Although this is somewhat smaller than the Mach num- 
ber of « 0.5 suggested by the resonance line studies of 



Perseus by Churazov et al. (2005), again our model only 
provides a lower limit on this value as it does not include 
all mechanisms that drive turbulence. 

Comparing the turbulent length scale between the two 
runs, we again find that the enhancement in the Sedov 
bubble case is not dramatic. Apart from a noticeable in- 
crease in L at R <; 20 kpc in the Sedov bubble runs, both 
runs achieve similar values of L, which averages » 5 kpc. 
Combining L and V to construct a diffusion coefficient 
gives typical values of 500 kpc km/s, which are in excel- 
lent agreement with the effective diffusion coefficient in- 
ferred by the abundance profiles studies of Perseus by Re- 
busco et al. (2005), again suggesting that more complete 
self-regulating models with subgrid turbulence could do 
well at explaining this and other clusters. 

In all cases, the total turbulent energy, shown in the 
the bottom panels, is « 1-2% of the buoyant energy in 
the bubbles, and even a smaller fraction of -©buoyancy + 
-©expand- Thus, as suggested above, turbulent dissipation 
is not likely to be an important factor. Rather it is the 
turbulent mixing of heated and enriched gas that plays a 
fundamental role in the evolution of AGN-heated galaxy 
clusters. 

6. CONCLUSIONS 

A wide range of observations suggest that AGN- 
feedback plays a key role in the evolution of cool-core 
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Fig. 18. — Radial profile of mean-flow quantities in Sedov 
multiple-bubble runs at times of 100 Myrs, 200 Myrs, and 300 
Myrs, arranged in columns from left to right. From top to bottom 
the rows show the average temperature, entropy, and metal 
density. Symbols are as in Figure 16. 

galaxy clusters. At the same time, theoretical studies 
have pointed out some of the many physical mechanisms 
that may be important in this evolution, including vis- 
cosity, magnetohydrodynamic effects, heat conduction, 
and cosmic-ray heating. While any or all of these may 
operate in nature, none of them can be fully understood 
without first accurately capturing the underlying hydro- 
dynamic evolution of the ICM. It is with this basic goal 
in mind that we have used a subgrid turbulence model to 
study AGN heating in an inviscid fluid, neglecting other 
effects such as magnetic fields and heat-conduction. 

Within this context, our study has been focused on 
two key issues: the growth of instabilities and turbulence 
caused by AGN-heated bubbles and the role of these in- 
stabilities in determining the evolution of the bubbles 
and the surrounding medium. Clearly, there are other 
sources of turbulence that might add similar levels of 
turbulent energy into the ICM, such as mergers of large 
subclusters or motions of galaxies within a cluster. Like- 
wise, the DT06 subgrid turbulence model that we have 
employed does not include all processes that generate 
turbulence, such as the shear-driven Kelvin-Helmholtz 
instability. However, our tests show that it is effective 
in capturing the growth of the extremely important RT 
and RM instabilities. Thus although many aspects of 
cluster evolution remain uncertain and beyond the scope 
of this work, there are a number of robust conclusions 
that we can make about the role of these two instabili- 
ties in shaping the evolution and impact of AGN-heated 
cavities in clusters. In particular we find that: 

• Many of the RT and RM unstable modes that 
drive the evolution of the bubbles evolve on scales 
that are far below the resolution limits of cur- 
rent simulations. The superposition of these un- 



stable modes smears out the interface between the 
bubbles and the ambient medium, transforming 
them into clouds of mixed material that stay in- 
tact and expand as they rise in the stratified clus- 
ter medium. This mixing can explain the coherent 
X-ray cavities detected in clusters of galaxies. The 
subgrid-turbulence model also greatly reduces the 
sensitivity of our results on the resolution of the 
computational grid and the detailed choice of ini- 
tial conditions. 

• Within the clouds, turbulent motions quickly at- 
tain typical velocities of ps 200 km s _1 , roughly 
10% of the internal sound speed and 20% of the 
sound speed of the surrounding ICM. Similarly, the 
scale of the turbulent eddies rises swiftly but does 
not exceed the ss 30 kpc scale of clouds. A typical 
turbulent diffusivity is then pa 500 km/s kpc, which 
is in excellent agreement with the diffusion coeffi- 
cient inferred by the abundance profiles studies of 
Perseus by Rebusco et al. (2005). 

• Subgrid turbulence is likely to enhance metal 
transport significantly. In our fiducial single- 
bubble evacuated pure-hydro run, metal transport 
is halted at ps 50 kpc by bubble disruption, which 
occurs when RT instabilities have shredded the 
evacuated region into resolution-limited cavities. 
Yet, in the subgrid-turbulence run, small-scale fluc- 
tuations act to keep the cloud coherent for much 
longer, thus causing it to distribute metals out to 
larger radii. 

• In the case where the bubbles produce shock 
waves, the turbulent velocities and length scales 
are systematically enhanced with respect to the 
evacuated-bubble run, due to the RM instability. 
While such ps 300 km/s velocities are roughly 50% 
greater than in the runs with evacuated bubbles, 
they are still well below the ps 500 km/s radial 
velocities of the clouds. However they are large 
enough to be probed by studying the emission lines 
of heavy ions with the future Constellation- X satel- 
lite, which will have an envisaged spectral resolu- 
tion of 1-2 eV. 

• Calculating the turbulent kinetic energy produced 
by the rising bubbles, we find that this is only pa 1% 
of the total energy available for the bubbles to heat 
the cluster. Hence, we do not expect the energy of 
the turbulent motions themselves to play a major 
role in the heating of the cool-core regions of the 
ICM. Rather, the main role of turbulence is to in- 
crease the efficiency with which the thermal energy 
of the rising clouds is mixed into the surrounding 
gas. 

• Finally, runs that include radiative cooling and 
multiple episodes of AGN-feedback indicate that 
the impact of turbulence continues to increase with 
successive generations of heating. This is because 
turbulence driven by previous feedback events re- 
mains behind, enhancing the mixing of subsequent 
bubbles. While in our pure-hydro runs, the bub- 
bles deposit most of their internal energy at the 
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Fig. 19. — Radial profile of turbulence quantities in the evacuated multiple bubble run (5DCER, solid lines), and the Sedov multiple- 
bubble run (5DCSR, dashed lines), at times of 100 Myrs, 200 Myrs, 300 Myrs, and 450 Myrs, arranged in columns from left to right. 
From top to bottom the rows show the average t urbu lent velocity and txurbulent length scale as a function of radius, and the total turbu- 
lent kinetic energy within this radius. As in Figure lMl the shaded regions in the lower panels show 1 — 2 of Sbuoyancy an d S cxpan( j+.E'| : > uoyancy . 



resolution-dependent radius at which they are dis- 
rupted, the turbulent motions captured by the sub- 
grid model lead to more gradual heating, corre- 
spondingly shallower temperature and entropy pro- 
files, and shallower metal gradients. 

In summary, the properties, evolution, and appearance 
of AGN-blown bubbles in clusters are substantially dif- 
ferent when one properly accounts for turbulent motions 
on scales well below the limits of current pure-hydro sim- 
ulations. Although accounting for these motions may not 
provide the ultimate solution to many of the mysteries 
surrounding galaxy clusters, it is a crucial step forward 
in the modeling of feedback and the interaction between 
AGN and the ICM. Subgrid models such as the one de- 
veloped in DT06 provide us with tools for capturing this 
physics. In fact, the numerical methodology presented 
here is likely to have applications in other areas of astro- 
physics where hydrodynamic modeling of RT instabilities 
is crucial, such as supernovae, supernova remnants, and 
galactic winds. It is clear that while many physical pro- 
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